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ABSTRACT 

\ Combustion of liquid fuels in the form of spray droplets is simulated numerically in this paper. 
Various, vaporization models are examined as to their performance in finite element calculations involving ~ 
turbulent, flow field. The Eulerianjdoordinate for the gas and Lagrangian Coordinate for the liquid spray ^ 
droplets af^coupled through source terms being updated in the equations of continuity, momentum, and 
energy. The f^and modified eddy breakup models are used for simulating turbulent spray combustion 
flow field. Numerical results for the droplet trajectories, droplet heating, recirculation characteristics, and 
effects of evapolation models are evaluated. It is also shown that the finite element method is 
advantageous in dealing with complex geometries, complex boundary conditions, adaptive unstructured 
grids. 

V\ 

1. INTRODUCTION 

One of the most critical aspects in spray combustion is an adequate vaporization model in which 
the effect of neighboring droplets on the rates of heat and mass transport and vaporization for any given 
droplet can be properly taken into account [1,2]. This is particularly important if the ambient gas 
temperature is high so that droplet life times and droplet heating times are of the same order of 
magnitude. Difficulties arise also in turbulent flow where the k — € model found satisfactory in 
nonreacting gas media may not be applicable in reacting spray combustion. A similar question is raised 
as to the adequacy of modified eddy breakup model to determine the mixing rates of reactants. The 
main issue in this paper, however, is not the development of such models, but rather the implementation 
of computational techniques for the turbulent spray combustion using the currently available physical 
models. Three vaporization models [2,3] are examined. 

We present finite element calculations [5] of turbulent spray combustion [6] using the Eulerian 
Coordinate for the gas and Lagrangian Coordinate for the liquid spray droplets, the k — € turbulence 
model, the eddy breakup model for predicting the mixing rates of the reactants. In turbulent reacting 
flows in the Eulerian coordinates, the velocity and pressure fields are strongly coupled with various source 
terms including turbulence, gas/liquid— phase interaction, and chemical reaction rates. 
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To prevent the spurious pressure modes, the mixed interpolation finite element method is used. Such a 
coupled solution eliminates the need for transformation of the continuity equation into a pressure or 
pressure correction equation in the sequential method. The coupled method used here is relatively 
insensitive to Reynolds numbers and grid aspect ratio. Other transport equations are solved sequentially. 
The solution procedure for gas— phase equations is similar to the single— phase turbulent reacting flows. 
The ordinary differential equations governing the droplet field in the Lagrangian coordinates are 
integrated using an explicit second— order Runge— Kutta method. The gas-phase properties at the 
characteristic location are calculated by linear interpolation of the four isoparametric finite element 
Eulerian nodal values in the computational element containing the droplet. The characteristic source 
terms at the Eulerian grid is evaluated by superimposing the nonlinear source term of each characteristic 
to the four surrounding grid points using the volume— weighted linear interpolation. 

Applications include calculations of droplet trajectories as well as the complete flow field variables 
in a centerbody combustor. It is shown that the present finite element model predicts the qualitative 
features of the turbulent spray combustion satisfactorily! pending verification by experimental 
measurements. The computational results show that the variations of thermophysical properties and the 
droplet heating model have the significant effects on the droplet history and the gas— phase flowfield near 
the injection region. 

2. GAS FIELD EQUATIONS IN EULERIAN COORDINATES 


The governing equations in Eulerian coordinates include equations of continuity, momentum, mass 
fractions, turbulent kinetic energy, dissipation, and concentration fluctuations. 
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with n^ — no. of droplets, and m^ = droplet vaporization rate. Here all variables are time — averaged 


mean quantities. 
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Turbulent Kinetic Energy 
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with (J ^ representing Schmidt or Prandtl number for the dependent variable <f). Note also that the 

gas— phase equations are non— dimensionalized using the inflow conditions of the combustor under study, 

leading to Re , the characteristic Reynolds number, M , the characteristic Mach number, and L, , the 
c c k 

ratio of the gas— phase length scale to the initial droplet radius of the k th group. 

The reaction rate R^ is determined from either the mixing rates of the reactants or the chemical 
reaction, whichever slower. The mixing rates of the reactants are obtained by the eddy breakup reaction 
model [7] assuming that the gas is composed of alternating fragments of unburned fuel and almost fully 
burned gas. The eddy breakup reaction rate is assumed to depend on the rate i/g at which the fragments 
of unburned fuel are broken into smaller fragments by the action of turbulence, and is taken to be 
proportional to the decay of turbulence energy e/k. 

The reaction rate may also be controlled by chemical kinetics when the mixing to the reactants is rapid 
[8]. Thus the reaction rate is expressed as 


R. = min 
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The specific heat of the mixture is given by 



Here by assuming equal binary diffusion coefficients for all the species and by knowing the mass fractions 
of fuel and oxidizer and the stoichiometric relationship of the given hydrocarbon— air mixture, the mass 


fractions of the remaining species (O , N , CO , and H 0) can be determined. The variation of the 

2 2 2 2 

specific heat, with temperature may be written as [9, 10]. 
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The density distribution in the gas phase is calculated from 


P = p/rt S Y./W. (12) 

The model constants used in the above equations are as follows: C = 1.44 , C = 1.92 li = 0.09 

1 2 

<7, = 1.0 , a - 1.217 , a = 0.9 , a = 0.9 , a, = 0.9 , C = 2.8 , C = 2.0 . 
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3. LIQUID SPRAY DROPLET FIELD EQUATIONS IN LA GRAN GIAN COORDINATES 


The liquid— phase equations are based on the Lagrangian formulation of the droplet trajectory, transient 
heating, and vaporization. The effective conductivity model used in the present study assumes the 
quasi— steadiness in the gas film. The instantaneous vaporization rate of the droplet is controlled by the 
transient process of heating of the liquid inside the droplet. In the limiting cases of the droplet 
impingement on the chamber walls the droplets undergo instantaneous vaporization. The governing 
equations for the liquid— phase are as follows: 
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Reynolds Number 
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Various evaporation models have been proposed and they 


applications. The droplet evaporation rate and the heat balance 
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appear to be quite sensitive in numerical 
equation may be expressed as 
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where m. is the evaporation rate; N and N denote the evaporation coefficient and the correction factor 

K Ad 

for convection effects, respectively; T, is the droplet surface temperature; Q is the heat transferred into 

K u 

the droplet interior; is the droplet mass; Cp^ is the specific heat of the droplet. We consider the 
following two models: 

(a) Model I [3] 


In this model the evaporation coefficient and the correction factor for convection effects are given by 
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where B is the mass transfer number. 


The heat energy, Q , takes the form 
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Here Y^ and are the mass fraction and the fuel vapor pressure at the droplet surface, respectively; P f 
the ambient pressure, is the sum of the fuel vapor pressure and the air partial pressure at the droplet 
surface; and W^. are the molecular weights of fuel and air, respectively. For any given value of 
surface temperature, the vapor pressure is readily estimated from modified Clausius— Clapeyron equation. 
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Values of a , a , and a are found in [3,4]. If the heat transfer to the droplet interior is neglected, then 

12 3 

we set ^ = H^. Also, to keep the heat transfer rate of vaporization (m H positive, it is 


v,eff 

necessary to maintain Y^ — 0. 


(b) Model 0 [3] 

Model II is the same as Model I except for the effective latent heat of vaporization. Considering the heat 

transferred in the droplet interior (Q^)» the effective latent heat of vaporization can be obtained from 
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where B is the heat transfer number 
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(c) 


Model m [4] 


This model includes the effects of variable thermophysical properties, non — unitary Lewis number in the 
gas film, the effect of the Stefan flow on heat and mass transfer between the droplet and the gas, and the 
effect of internal circulation and transient liquid heating. Thus, the evaporation coefficient and the 
correction factor for convection effects are as follows: 
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The heat energy takes the form 
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The average properties of air— vapor mixture may be determined at the following reference 
temperatures and compositions [11]: 
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In the limiting cases of the droplet impingement on the chamber walls, the droplets undergo 
instantaneous vaporization. A similar treatment of the droplet is considered when approximately 97% of 
the mass of the droplets is vaporized. In case of the droplet passage through the plane of symmetry, 
another droplet with similar instantaneous properties and physical dimensions, but with the mirror image 
velocity vector is injected into the flow field, 

4. COUPLING OP LIQUID AND GAS PHASES 

The liquid— phase equations for spray droplets are advanced in time by an explicit second— order 
Runge — Katta method using time steps much smaller than the gasr— phase equations. Based on the known 



gas— phase properties, the liquid— phase equations are first advanced in time from the nth time level to 
the (n-hl)th time level corresponding to the gas— phase time step Atg as follows: 

(1) Interpolate linearly the gas— phase properties from the fixed Eulerian 

locations to the characteristic location. 

(2) Integrate the liquid— phase equations over a time step, At 

(3) Redistribute the source terms evaluated at the characteristic location among the 
Eulerian nodes surrounding the characteristic. 

Steps (1) — (3) are repeated until the liquid— phase equations are advanced over a time 
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step, Atg. 

Solve the gas— phase equations using the time implicit scheme. 

Repeat steps (1) — (5) until the iteration converges before advancing to the next step. 


The injected spray and the droplet flow is assumed to be conical such that 


u, = u. cos 0 , 
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v, = v, sin 6 
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where 6 is the half— cone angle. The mass flow rate of the fuel is determined from the stoichiometric 

conditions and the injection time interval for a new characteristic to appear is based on the Eulerian 

mesh size Ax and injection velocity of the droplet. 
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Thus the number of droplets in each characteristic is given by 
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Once the droplets are injected its subsequent behavior is determined by the governing equations. 
5. BOUNDARY CONDITIONS 


Steep gradients and a relatively low level of turbulence prevail along the wall. To account for rapid 
changes in velocity profile close to the wall, the so-called wall function method is employed [12], In the 
context of finite elements, we assume that the shear stress is constant up to a distance 5 from the wall 



such that 


r 

w 


with 


u 

* 7 

p ku C« 

u 

&iE<5 + 


for 8 <11.6 
for 8* > 11.6 


♦ pS cj k* 

s = 

Once the near— wall values of the shear stresses are evaluated, we calculate the near— wall values of 
the k and e as follows: 
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The near— wall heat flux is determined by 
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where the function P(pr) is of the form. 
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ff ere r w anc ^ are specified as Neumann boundary conditions in the momentum and energy equations. 

6. APPLICATIONS 


Numerical results include comparisons of three different models of vaporization models applied to a 
typical turbulent spray combustion. Figures l.a through l.d show the vaporization characteristics of 
n-decane droplets of initial radius r^ = 50 fJm and T k q = 300°K which are injected into the air 
stream at P x = lOatm, = 1500°K, and Au q = u - U k = lOm/s. The temporal variation of the 
non-dimensional droplet surface area, surface temperature, vaporization rate (m/m ), effective heat 
transfer rate of evaporation (mH^ ,efP' are shown in Figures l.a, l.b, l.c, and l.d, respectively. Here m 
is the initial droplet mass. 



The numerical results indicate that the droplet heating process takes a considerable part of the droplet 
lifetime. Note here that the vaporization histories of Model 1 are the same as those of Model 2 except 
when the effective vaporization heat transfer rate is used for Model 2. The transient heating period of 
Model 1 or Model 2 is much shorter than Model 3. At the final stage, the droplet surface temperature 
for both models approaches an equilibrium value which is the wet— bulb temperature. Model 1 or 
Model 2 leads to a higher wet-bulb temperature than Model 3. Since the transfer number of Model 3 is 
based on the effective latent heat vaporization, Model 3 has a slower vaporization than Model 1 or 
Model 2 in the initial heating period. Figure l.d shows the temporal variations in the effective heat 
transfer vaporization (mH ,efP‘ Models 2 and 3 have the large heat transfer rate in the early 
evaporation period which is controlled by the effective latent heat vaporization. However, the heat 
transfer rate of vaporization for Model 1 is controlled by the vaporization rate because the effective latent 
heat in Model 1 is obtained by neglecting the heat transfer to the droplet interior. Since Model 2 has a 
faster vaporization than Model 3 in the droplet heating period, Model 2 shows a larger effective heat 
transfer rate than Model 3. 

Figures 2. a through 2.d show the vaporization characteristics of n— decane droplets in the air— fuel 
mixture medium with Y foo = 0.4. All other conditions are the same as the previous example. This 
situation is frequently encountered in the actual spray combusting flows. The vaporization history of 
Model 1 is the same as the previous case because this model does not account for the effect of the 
ambient fuel vapor. The introduction of any cold droplet into a surrounding of its own vapor results in 
vapor condensation on the droplet surface. In the condensation period, the mass transfer number, B^, 
becomes negative if the mass fraction of fuel at the interface is less than the mass fraction of fuel at the 
ambient. As a result, the droplet undergoes the increase in radius. It can be seen that the droplets with 
its ambient fuel vapor evaporates much faster than the previous case without the ambient fuel vapor. In 
fuel vapor environment, the droplet Reynolds number increases due to variable property effects related to 
the liquid dynamic viscosity. The effect of an increase in Reynolds number on heat and mass transfer 
obviously causes the Nusselt number and Sherwood number to increase. Therefore, the higher droplet 



Reynolds number results in increase of the vaporization rate. During the condensation period, Model 2 
yields a larger condensation rate than Model 3. Shortly after the condensation period, Model 2 exhi bit! a 
faster vaporization than Model 3. However, at the final stage of vaporization, Model 2 shows a slower 
evaporation than Model 3. Model 2 yields a larger effective heat transfer rate than Model 3 during the 
entire vaporization period except for the condensation period and the final stage. 

The next example is for the computation of the spray combusting flows using three vaporization 
models with the geometry of a centerbody combustor shown in Figure 3 and the initial and boundary 


conditions summarized in Table 1. The dimensions of the combustor are the same as those used in Raju 

and Sirignano [13]. However, the present study uses the variable thermophysical properties. The 

injected spray is assumed to comprise four conical streams and half-angles of the corresponding streams 

at ^ = 5, 15, 25, and 35 degrees. The computations are performed for three vaporization models 

discussed earlier. In the limiting cases of the droplet impingement on the chamber walls, the droplets 

undergo instantaneous vaporization. Similar treatment of the droplet is considered when 97% of the 

mass of the droplet is vaporized. In case of the droplet passage through the plane of symmetry, another 

droplet with similar instantaneous properties and physical dimensions but with the mirror image velocity 

vector is injected into the flowfield. The time steps for the steady state are:At. . = 1.6 m/s At =16 

inj ’ g 


ms, and At, = 0.04ms. 
l,m 

Figures 4.a through 4.c show droplet trajectories and vaporization processes for three vaporization 
models. The four droplet groups can be identified by the volume of the droplet and the characteristic 
location. It is seen that the droplet motion is initially governed by the droplet inertia force before the 
drag force causes the droplets to decelerate and the droplet path is eventually determined by the 
gas— phase flow field. Most of the vaporization occurs within the recirculation zone because the smaller 
droplets are unable to penetrate downstream. Because of the strong negative radial gas— phase velocity 
field near the injector, the droplet trajectories are significantly affected by the gas— phase velocity field 
especially for the droplet characteristic with the lowest injection angle, 0 = 5 degrees. 



The strong negative radial gas— phase velocity field in the injector region results from the large drag force 
term interacting with the source terms of radial momentum equation. Model 1 has a faster evaporation 
rate than two other models. Slight differences between Model 2 and Model 3 exist in the droplet 
trajectories corresponding to the final stage of vaporization. 

Velocity vectors for three vaporization models are shown in 5.a through 5.c. The velocity field for 
three models has the similar secondary recirculation zone. This is due to the gas— droplet interaction in 
recirculation zone having the high vaporization rate. Slight differences among three models exists in the 
downstream side of the recirculation region. 

Figure 6 shows contours of temperature for three vaporization models. The temperature difference 
between two adjacent lines is about 150°K. The maximum and minimum temperatures of the gas-field 
are about 2800°K and 700°K. The low temperature near the injector results from the cooling effect of 
the vaporization process. This region is also characterized by large temperature gradients. 

Radial profiles of temperature for three vaporization models are presented in Figure 7. Since the 
effective latent heat for Model 1 predicts a much higher temperature than the other two models in the 
fuel— rich side of three axial locations, Model 1 exhibits the highest temperature, followed by Model 3 and 
Model 2. Temperature distributions in the fuel — rich region mainly depend on the effective heat transfer 
rate of vaporization (m ^). 

Figure 8 shows contours of the fuel mass fraction for three vaporization models. The values of the 
contour lines vary from 0.005 to 0.655 with the mass fraction difference of 0.05. The large concentration 
of fuel vapor in the recirculation region is due to insufficient mixing of fuel and air. 

Finally, radial profiles of the fuel mass fraction for three vaporization models are shown in 
Figure 9. Since Model 1 predicts a faster evaporation due to a higher temperature near the injector 
region (x = 0.004m), Model 1 yields a much larger fuel mass fraction than the other two models. At 
some distance from the injector (x = 0.08m, 0.12m), Model 1 predicts the lowest mass fuel fraction 
resulting from the shortest droplet lifetime with the high gas— temperature environment. Model 2 shows 
a larger fuel mass fraction than Model 3 at some distance downstream (x = 0.08m, 0.12m), because 
Model 2 yields a faster evaporation related to the complex effects of variable properties. 


7. 


CONCLUSIONS 


Numerical analysis using the finite element Eulerian— Lagrangian approach for various vaporization 
models for turbulent spray combustion has been shown to be effective. The following comments are 
offered: 

(1) Variations of the thermophysical properties and the droplet vaporization models are sensitive to the 
droplet histories and the gas— phase flowfield especially close to the fuel injector region. 

(2) Droplet trajectories are greatly influenced by the choice of existing vaporization models which 
determine the vaporization rate and the effective latent heat of vaporization. 

(3) Model 1 has a faster vaporization rate, leading to the droplet trajectories affected more rapidly by the 
gas— phase flowfield than Models 2 and 3. 

(4) The gas— phase velocity field for all vaporization models appears to have a similar secondary 
recirculation zone. The gas— droplet interactions play a negligible role for the formation of secondary 
recirculation zone. 

(5) Temperature distributions near the fuel— rich injector are significantly affected by the transient 
droplet heating in terms of the effective latent heat of vaporization. Model 1 predicts the highest 
temperature distributions, followed by Model 3 and Model 2. 

(6) Close to the fuel injector, Model 1 yields a much larger fuel mass fraction than the other two models, 
reversing the trend downstream because of the shortest droplet lifetime with the higher gas— phase 
temperature distribution. 

(7) Model 2 predicts a slightly larger fuel mass fraction than Model 3 downstream, due to a faster 
evaporation rate. 

(8) In terms of computational strategies, the finite element method would be convenient in dealing with 
complex geometries, boundary conditions, and adapative unstructured grids. 
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Fig 1 Vaporization Characteristics on n— decane droplets injected into air flow. 
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FIG 2 Vaporization Characteristics of n-decane droplets in the Air-fuel Mixture 
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FIG 6 Contour* of Temperature for Three Vaporisation Model*, 

Contour lateral 50°K, Minimum 700°K, Maximum 2800°K 


FIG 7 Radial Profile* of Temperature for Three Vaporisation Modal* 
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FIG 8 Fuei Mail Fraction*, Contour Intern! 0.05, Minimum 0.05, Maximum 0.855 


FIG 9 Fuel Maw Fractions in Radial Direction 






